#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include <string.h>
#include <unistd.h>

#define SAMPLE_COUNT 200    // 曲线采样点数
#define TAU_COUNT 7         // 张力τ的数量
#define MATRIX_SIZE 4       // Catmull-Rom基矩阵维度

// ========================= 全局常量 =========================
// 控制点 (x,y)，与Python一致
static const double P[MATRIX_SIZE][2] = {{1.0, 0.0}, {2.0, 3.0}, {3.0, -1.0}, {4.0, 5.0}};
// 张力τ取值，与Python一致
static const double taus[TAU_COUNT] = {-1.5, -1.0, -0.5, 0.0, 0.5, 1.0, 1.5};
// Gnuplot颜色（对应viridis配色）
static const char* colors[TAU_COUNT] = {
    "#440154", "#482878", "#3E4A89", "#31688E",
    "#26828E", "#1F9E89", "#35B779"
};

// ========================= 矩阵/向量运算工具 =========================
/**
 * 构建Catmull-Rom基矩阵 M(τ)
 * @param tau 张力参数
 * @param M 输出：4x4基矩阵
 */
void catmull_rom_matrix(double tau, double M[MATRIX_SIZE][MATRIX_SIZE]) {
    // 按Python公式填充矩阵
    M[0][0] = 0.0;        M[0][1] = 1.0;        M[0][2] = 0.0;         M[0][3] = 0.0;
    M[1][0] = -tau;       M[1][1] = 0.0;        M[1][2] = tau;         M[1][3] = 0.0;
    M[2][0] = 2*tau;      M[2][1] = tau - 3.0;  M[2][2] = 3.0 - 2*tau; M[2][3] = -tau;
    M[3][0] = -tau;       M[3][1] = 2.0 - tau;  M[3][2] = tau - 2.0;   M[3][3] = tau;
}

/**
 * 向量-矩阵乘法 (1x4) · (4x4) = (1x4)
 * @param vec 输入向量（长度4）
 * @param mat 输入矩阵（4x4）
 * @param res 输出结果（长度4）
 */
void vec_mat_mult(const double* vec, const double mat[MATRIX_SIZE][MATRIX_SIZE], double* res) {
    memset(res, 0, MATRIX_SIZE * sizeof(double));
    for (int i = 0; i < MATRIX_SIZE; i++) {
        for (int j = 0; j < MATRIX_SIZE; j++) {
            res[i] += vec[j] * mat[j][i];
        }
    }
}

/**
 * 向量-点矩阵乘法 (1x4) · (4x2) = (1x2)
 * @param vec 输入向量（长度4）
 * @param pts 输入点矩阵（4x2）
 * @param res 输出结果（x,y）
 */
void vec_pts_mult(const double* vec, const double pts[MATRIX_SIZE][2], double* res) {
    res[0] = 0.0; res[1] = 0.0;
    for (int i = 0; i < MATRIX_SIZE; i++) {
        res[0] += vec[i] * pts[i][0];
        res[1] += vec[i] * pts[i][1];
    }
}

// ========================= Catmull-Rom曲线计算 =========================
/**
 * 计算指定τ的Catmull-Rom曲线点
 * @param tau 张力参数
 * @param u_vals 参数u向量（0~1，长度SAMPLE_COUNT）
 * @param x_out 输出x坐标数组（长度SAMPLE_COUNT）
 * @param y_out 输出y坐标数组（长度SAMPLE_COUNT）
 */
void catmull_rom_curve(double tau, const double* u_vals,
                       double* x_out, double* y_out) {
    // 1. 构建基矩阵M
    double M[MATRIX_SIZE][MATRIX_SIZE];
    catmull_rom_matrix(tau, M);

    // 2. 遍历所有u值计算曲线点
    for (int i = 0; i < SAMPLE_COUNT; i++) {
        double u = u_vals[i];
        // 构建参数向量U = [1, u, u², u³]
        double U[MATRIX_SIZE] = {1.0, u, u*u, u*u*u};
        // 计算C = U·M
        double C[MATRIX_SIZE];
        vec_mat_mult(U, M, C);
        // 计算曲线点 = C·P
        double pt[2];
        vec_pts_mult(C, P, pt);
        // 保存结果
        x_out[i] = pt[0];
        y_out[i] = pt[1];
    }
}

// ========================= Gnuplot绘图 =========================
/**
 * 初始化Gnuplot管道
 * @return Gnuplot文件指针
 */
FILE* gnuplot_init() {
    FILE* gp = popen("gnuplot -persist", "w");
    if (!gp) {
        perror("无法启动gnuplot");
        exit(EXIT_FAILURE);
    }
    // 设置绘图参数
    fprintf(gp, "set title 'Catmull‑Rom 样条（不同张力 τ）'\n");
    fprintf(gp, "set xlabel 'x'\n");
    fprintf(gp, "set ylabel 'y'\n");
    fprintf(gp, "set grid lw 1 lt '--' alpha 0.6\n");
    fprintf(gp, "set key right bottom\n");
    fprintf(gp, "set terminal wxt size 800,500\n");
    return gp;
}

/**
 * 绘制控制点
 * @param gp Gnuplot文件指针
 */
void plot_control_points(FILE* gp) {
    fprintf(gp, "plot '-' with points pt 7 ps 1.5 lc rgb 'black' title '控制点', \\ \n");
    for (int i = 0; i < MATRIX_SIZE; i++) {
        fprintf(gp, "%.6f %.6f\n", P[i][0], P[i][1]);
    }
    fprintf(gp, "e\n");
    fflush(gp);
}

/**
 * 绘制指定τ的曲线
 * @param gp Gnuplot文件指针
 * @param tau 张力参数
 * @param color 曲线颜色
 * @param x 曲线x坐标数组
 * @param y 曲线y坐标数组
 * @param idx 曲线索引（用于拼接绘图命令）
 */
void plot_curve(FILE* gp, double tau, const char* color,
                const double* x, const double* y, int idx) {
    char buf[1024];
    // 拼接绘图命令
    if (idx == 0) {
        snprintf(buf, sizeof(buf),
                 "'-' with lines lw 2 lc rgb '%s' title 'tau=%.1f'",
                 color, tau);
    } else {
        snprintf(buf, sizeof(buf),
                 "'-' with lines lw 2 lc rgb '%s' title 'tau=%.1f' \\ \n",
                 color, tau);
    }
    fprintf(gp, "%s", buf);

    // 写入曲线数据
    for (int i = 0; i < SAMPLE_COUNT; i++) {
        fprintf(gp, "%.6f %.6f\n", x[i], y[i]);
    }
    fprintf(gp, "e");
    if (idx < TAU_COUNT - 1) {
        fprintf(gp, ", \\ \n");
    } else {
        fprintf(gp, "\n");
    }
    fflush(gp);
}

// ========================= 主函数 =========================
int main() {
    // 1. 生成参数u向量（0~1，200个采样点）
    double u[SAMPLE_COUNT];
    double step = 1.0 / (SAMPLE_COUNT - 1);
    for (int i = 0; i < SAMPLE_COUNT; i++) {
        u[i] = i * step;
    }

    // 2. 初始化Gnuplot
    FILE* gp = gnuplot_init();

    // 3. 绘制控制点（先绘制基础图层）
    plot_control_points(gp);

    // 4. 计算并绘制所有τ对应的曲线
    for (int i = 0; i < TAU_COUNT; i++) {
        double tau = taus[i];
        const char* color = colors[i];

        // 计算曲线点
        double x[SAMPLE_COUNT], y[SAMPLE_COUNT];
        catmull_rom_curve(tau, u, x, y);

        // 绘制曲线
        plot_curve(gp, tau, color, x, y, i);
    }

    // 5. 关闭Gnuplot管道
    pclose(gp);

    printf("绘图完成！请查看弹出的Gnuplot窗口。\n");
    return 0;
}